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Abstract 

Recently, we introduced a new test for distinguishing regular from chaotic 
dynamics in deterministic dynamical systems and argued that the test had cer- 
tain advantages over the traditional test for chaos using the maximal Lyapunov 
exponent. 

In this paper, we investigate the capability of the test to cope with moderate 
amounts of noisy data. Comparisons are made between an improved version 
of our test and both the "tangent space" and "direct method" for comput- 
ing the maximal Lyapunov exponent. The evidence of numerical experiments, 
ranging from the logistic map to an eight-dimensional Lorenz system of differ- 
ential equations (the Lorenz 96 system), suggests that our method is superior 
to tangent space methods and that it compares very favourably with direct 
methods. 
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1 Introduction 

Given time series data from a deterministic dynamical system, it is often of interest 
to determine whether the underlying dynamics are regular or chaotic. For example, 
heart rate data measured in electrocardiograms are believed to be chaotic for the 
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case of healthy patients but show regularity in the case of congestive heart failure 
(see Wagner & Persson jTU] for an overview). 

The usual test of whether a deterministic dynamical system is chaotic or non- 
chaotic is the calculation of the maximal Lyapunov exponent A. Standard references 
include [TJ HJ UJ |Hl EU UJJ • A positive maximal Lyapunov exponent indicates chaos: 
if A > 0, then nearby trajectories separate exponentially and if A < 0, then nearby 
trajectories stay close to each other. This approach has been widely used for dy- 
namical systems whose equations are known. If the equations are not known or one 
wishes to examine experimental data, then A may be estimated using the phase space 
reconstruction method of Takens ^B] or by approximating the linearisation of the 
evolution operator. 

In a previous paper [H], we introduced a binary test for distinguishing between 
regular and chaotic dynamics in deterministic dynamical systems. The test has two 
main advantages over computing the maximal Lyapunov exponent: 

(i) The test applies directly to time series data, so that phase space reconstruction 
is not required. Moreover, the form and nature of the underlying dynamical 
system is irrelevant; the test applies equally well to continuous time systems and 
discrete time systems, to experimental data and maps, to ordinary differential 
equation and partial differential equations. 

(ii) It is a binary test (in principle, the test yields or 1), so a numerically computed 
value of 0.01 say yields a definite conclusion, whereas such a value for the 
maximal Lyapunov exponent would not. 

The theoretical basis for the validity of any numerical test for chaos relies on 
the availability of unlimited noiseless data. In practice, it is necessary to work with 
limited amounts of contaminated data. The only way to determine the utility of the 
test in such situations is to try it out on different examples and see how it fares in 
comparison with existing methods. Such a comparison is carried out in this paper. 

In Section |2J we describe an improved version of the test introduced in In 
Section |21 we briefly review the tangent space and direct methods for computing 
maximal Lyapunov exponents. In Section 01 we apply the test to the logistic map 
contaminated with measurement noise. Here, our test compares extremely favourably 
to the tangent space method for computing the Lyapunov exponent. In Section 
the test is applied to contaminated data from an eight-dimensional ODE and is seen 
to compare favourably even to direct methods. 
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2 The (modified) test for chaos 



Our test for chaos in this paper is a slightly modified version of the test that we 
introduced previously [5]. In this section, we describe the modified test and compare 
it with the one in We focus throughout on discrete data sets 0(n) sampled at 
times n = 1,2,3,... (The continuous time case is similar, cf. [5].) Here, <p(n) is a 
one-dimensional observable obtained from the underlying dynamics. 

Choose a constant celat random. In practice, we choose a number of different 
values of c as detailed below. For each value of c, define 

n 

p(n) = cos t? c )> n = 1,2,3,... (2.1) 

i=i 

Next, define the mean square displacement 

1 N 

M{n) = lim — V[p(j + n) - p(j)] 2 , n = 1, 2, 3, . . . 

N^oo iv ' J 

i=i 

If the dynamics is regular (periodic or quasiperiodic), then with probability one M(n) 
is a bounded function of n. However, if the dynamics is chaotic (in a fairly mild sense), 
then with probability one M(n) = Vn + 0(1) for some V > 0. (For justifications of 
these statements, see |3] and the references therein.) Define the asymptotic growth 
rate of the mean square displacement 

K = lim 

rwoo log 71 

Then = signifies regular dynamics whereas K = 1 signifies chaotic dynamics. 
Next suppose that we have a finite amount of data 4>(n), 1 < n < N. Define 

= [pO' + - pU)} 2 - 

3=1 

Provided n it is reasonable to expect that M(n) scales with n in a similar way to 
before. To avoid logarithms of negative numbers, we plot log(M(n) + 1) against logra 
for 1 < n < Ni for some choice of Ni, 1 < iVi C iV. In practice, we have found it 
convenient to choose Ni = N/10. This choice is made throughout the paper without 
further comment. (A disadvantage of our test is that we cannot use the full available 
data set for our statistics in the computation of the mean square displacement.) We 
define K to be the slope of the line of best fit on the resulting data set, using least 
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square regression. The test is now that K close to zero signifies regular dynamics and 
K close to 1 implies chaotic dynamics. 

For short data sets (with or without noise) there is an inconvenient resonance phe- 
nomenon, where a given choice of c may resonate with frequencies in the underlying 
dynamics. For example, if the signal is 27r-periodic and we choose c to be an integer, 
then a simple argument using the Fourier series for 4>(n) shows that typically p(n) will 
grow linearly yielding K = 2. A near- integer choice of c results in K = eventually 
but the convergence is very slow; for a small number of iterates the computed value 
of K is likely to be closer to 1 than to 0. The situation for quasiperiodic dynamics 
is even worse since "bad" choices of c are now dense. Nevertheless, there is zero 
probability of making a bad choice in theory. In practice, numerical experimentation 
shows that even for short data series, the bad choices of c are rare but do occur from 
time-to-time. 

Our resolution of this problem is to choose several values of c at random, comput- 
ing K for each choice of c. Then we take the median value of K. (We do not take 
the average of K, since when c does fail it can fail quite badly.) 

To demonstrate the improvement in the modified test, we consider the logistic 
map 



varying the parameter \i in the range 3.5 < /i < 4 in increments of 0.001. Starting 
with initial condition 0.0001, we use N = 1,000 iterates, after a transient of 20,000 
iterates. We take 0(n) = x n . 

In Figure Q we show how the old method [S] compares with the modified test, 
where for the modified test we take the median of 100 different values for c. It is 
clear that there is a dramatic improvement in our test for chaos. (A comparison 
with traditional Lyapunov exponent methods is carried out in Section HJ) An obvious 
criticism is that K is always far from 1 for such short data sets. In later sections, 
we demonstrate that this issue is of little consequence in comparison with traditional 
methods. Our purpose in this section is solely to compare the modified test with our 
old test. 
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(a) (b) 

Figure 1: Plots of K versus \x for the logistic map f(x) = fix(l — x), with 3.5 < \i < 4 
using 1,000 iterates and noise-free data, (a) The old test [S]. (b) The new test with 
100 values of c. 
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In Figure 121 we compare the effects of taking 1, 10, 100 and 1000 different values 
of c. It is clear from these results that taking one choice of c is not effective and that 
taking 10 values of c is a dramatic improvement, though some periodic windows are 
still poorly defined. This is remedied by going to 100 values of c, and increasing to 
1000 does not make a significant further improvement. In the remainder of the paper, 
we use 100 values of c. 




Figure 2: Plots of K versus fi for the logistic map f(x) = fix(l — x), with 3.5 < fi < 4 
using 1, 000 iterates and noise- free data. The median value of K is taken from m 
values of c where (a) m — 1, (b) m — 10, (c) m = 100, (d) m = 1000. 



Remark 2.1 There are two differences between the modified test presented here 
and the test introduced in [3]. The first is that we now use several values of c. The 
second is that in [Hj we defined p(n) by the more complicated prescription: 

0(n + l) =6(n)+c + ^(n) 1 > 22) 

p(n + 1) = p(n) + 4>(n) cos(#(n)) (' 

We note that the modified p(n) in equation (j2.1j) can be recovered by removing the 
4>(n) term in the formula for 9(n + 1) in equations (J2.2|) . 



6 



Numerical experimentation (specifically with the forced van der Pol oscillator 
considered in [5]) shows that the modified definition of p in (|2.1j) is inferior for long 
data sets with no noise when only one choice of c is used. This was our reason for 
including the extra 4>(n) term in (|2.2j) . However, this term is no longer advantageous 
when 100 choices of c are sampled. Moreover, (j2.1|) now has the advantage that p(n) 
scales linearly with the data and so is far less susceptible to measurement noise than 
the test in 
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3 Traditional methods 



Suppose that 4>{n) is the data set originating from an unknown deterministic dynam- 
ical system. (As usual <fi is some one-dimensional observable of the system under 
investigation.) To compute the maximal Lyapunov exponent, it is first necessary to 
reconstruct the dynamics [IB]. Define the m-dimensional delay vector 

a = {(j)(n), <f>{n + t), <f>{n + 2r), • • • , <f>{n + (m - l)r)}. 

Here, the delay time r > is chosen to be an integer. If the underlying dynamics x„ 
lies inside a <i-dimensional phase space, then the delay reconstruction map x n i— ► £ n is 
an embedding provided m > 2d+ 1 (corresponding to Whitney's embedding theorem 
PU]). Sauer et al. [TH H2] extended Takens' work by showing that it suffices to take 
m > 2d (A) where d$(A) is the fractal box-counting dimension of the attractor. Under 
these conditions, the reconstructed dynamics £i, • • • in M m faithfully represents the 
underlying dynamics. 

The maximal Lyapunov exponent for a map £ n+ i = /(£„) is defined to be 

A= hm Ilog(||D/g||). (3.1) 

There are (at least) two main numerical approaches for computing the maximal Lya- 
punov exponent: the tangent space method and the direct method. 

In tangent space methods, a model is fitted to the data to approximate the Ja- 
cobian. By (j3.1j) . the maximal Lyapunov exponent is defined as the product of the 
Jacobian of the the linearised flow along the given trajectory. If the underlying equa- 
tions are not known the Jacobian can be approximated. Such tangent space methods 
were first introduced by Eckmann & Ruelle and then further developed by Sano 
& Sawada ^3] and Eckmann et al. j3]. 

Direct methods use the fact that nearby trajectories separate on average asymp- 
totically at rate e Xn . The original algorithm by Wolf et al. allows the computation 
of the whole spectrum of Lyapunov exponents. Rosenstein et al. ^2] and Kantz jB] 
concentrated attention on the maximal Lyapunov exponent, and developed a more 
robust method. In the method by Rosenstein et al. ^2] for each m-dimensional vector 
£ in the reconstructed phase space, its nearest neighbour is determined. Rosenstein 
et al. calculate C(k) =< log \f k (£,)~ f h {£,*)\ > where the angles denote averaging over 
all £ = £i,£2> • • ■ The function C(k) shows roughly three different regimes. An initial 
regime of flat increase, a subsequent interval with exponential behaviour, and finally 
a plateau (because the separation cannot go beyond the extension of the attractor). 
The maximal Lyapunov exponent is determined by the slope of C(k) in the usually 
quite short range of exponential behaviour. 
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Remark 3.1 There are a number of inherent problems of both the tangent space 
method and the direct method, linked to phase space reconstruction. These difficul- 
ties, which include the choice of embedding dimension m and delay parameter r, are 
well-documented and we refer to Casdagli et al. j2] and Schreiber & Kantz ^H] for 
detailed discussions. 

4 Comparison with tangent space methods 

In this section, we compare our test with the tangent space method. Note that the 
definition of the maximal Lyapunov exponent (|3.1|) is itself a tangent space method. 
We use the method developed by Sano & Sawada [T3|. In particular we study the 
influence of measurement noise. 

As discussed in Section tangent space methods seek to approximate the Jaco- 
bian Df n (xo). In the case of measurement (additive) noise, the noise is amplified by 
higher order nonlinearities. Similarly, if the dimensionality of the underlying dynam- 
ical system is large, the evaluation of the diagonalised Jacobian requires extensive 
multiplication of the individual matrix elements of the Jacobian, and noise is again 
amplified by this process. 

Neither of these problems arises in our test. Measurement noise only enters our 
diagnostic variable p(n) linearly via the observation 4>(n) irrespective of the underlying 
dynamical system. 

We consider an illustrative example. Our contaminated data has the form <f>(n) = 
4>(n) + (A/yi00)?7 n where x n is the clean data, Af is the noise-level in percent, and 
771 , 772, .. . are i.i.d. random variables drawn from a uniform distribution on [—1,1]. 
We note that the results are similar in the case of normally distributed noise. 

As in Section |2l we study the logistic map 

x n+l = fix n (l - x n ), (4.1) 

varying the parameter fi in the range 3.5 < fi < 4 in increments of 0.001. We again 
use data sets consisting of iV = 1,000 iterates, after a transient of 20,000 iterates 
starting from the initial condition 0.0001. Throughout, we take 0( Tl^j — Xfi clS the 
observable. 

Our results for noise-free data are shown in the first column of Figure El As a 
benchmark, we compute the "exact" Lyapunov exponent making use of the explicitly 
given form of the map in (|4.1|) . see Figure Efa), Next, we use the tangent space 
method proposed by Sano & Sawada [12] for approximating the dynamics on the 
tangent space using phase-space reconstruction, see Figure Efb). The results from 
our modified test (Section^ with 100 values of c are shown in Figure HD^c). 
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Now we add measurement noise with noise-level M = 1%. As can be seen from 
Figure OJe) , the phase space reconstruction method by Sano & Sawada [12] performs 
rather poorly. Moreover, the numerical values and the overall qualitative behaviour of 
the maximal Lyapunov exponent with respect to varying \x depends very sensitively 
on the particular choice of the embedding dimension m and the value chosen to define 
nearest neighbours. In contrast, our test performs as well as it did in the noise-free 
case, see Figure Off)- 

It turns out that the tangent space method copes poorly with measurement noise 
even when phase space reconstruction is not required. To see that this is the case, 
we fed contaminated data into the exact expression for the Jacobian computed ana- 
lytically from (J4.1)) . In Figure 01 we show the results for the "exact" tangent space 
method compared with our test, both with 10% measurement noise. This seems to 
be conclusive evidence that our test is better than the tangent space method. We 
have also verified that our test deals comfortably with 20% noise. 

There is another advantage of our method, even in the noise-free case. Suppose 
that a whole scan through the bifurcation parameter is not available but instead data 
is only available for one particular value of the bifurcation parameter. A value of 
0.2 for the maximal Lyapunov exponent is inconclusive (see for example the periodic 
windows in Figure EJ), whereas a value of K = 0.01 indicates regular dynamics. In 
this sense, our test is a better absolute test than the Lyapunov exponent. 
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Figure 3: Plots of the Lyapunov exponent and K versus \x for the logistic map (|4.1|) 
with 3.5 < yU < 4 using 1, 000 iterates. First row: "Exact" Lyapunov exponent. 
Second row: Method by Sano & Sawada with m = 2 and 0.000001 as an allowed 
distance to define nearest neighbours. Third row: Our test with 100 different values 
of c. The first column uses noise- free data. The second column incorporates 1% 
measurement noise. 
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(a) (b) 



Figure 4: Plots of the Lyapunov exponent and K versus [i for the logistic map (|4.1|) 
with 3.5 < fi < 4 using 1,000 iterates with 10% noise, (a): "Exact" Lyapunov 
exponent, (b): Our test with 100 different values of c. 
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5 Comparison with direct methods 



In this section, we compare our test for chaos against the direct method for computing 
the maximal Lyapunov exponent. In particular, we use the direct method proposed 
by Rosenstein et al. |T2"] . 

The standard Lorenz equation is a 3-dimensional truncation of a certain system 
of partial differential equations (PDEs). This is too simple a system to present a 
challenge for testing for chaos. Hence, we study the n-dimensional Lorenz system, 
known as the Lorenz 96 system: 

^ = Xi_i(x i+ i - Xi_ 2 ) - Xi + r with z = l,---,n. (5.1) 

This system of ODEs was first introduced by Lorenz as an idealised model for the 
atmosphere 0, and has been studied by Orrell & Smith [TOj. In the following we 
choose n = 8. For 3.8 < r < 3.9 there are periodic windows, whereas for 5.4 < r < 
6.8 most regular windows are due to quasiperiodic dynamics [12]. Throughout, we 
integrate the system using a time step of 0.05 and calculate until 250, 000 units of 
time after a transient of 75, 000 units of time. However, we use only the data recorded 
after each 2.5 units of time, yielding a time series of N = 10, 000 data points. In this 
way, we obtain a discrete time series mimicking the situation for experimental data. 
As an observable we take <fi = x 2 + x 3 + x 4 , so 0(n) = x 2 (t) + x 3 (t) + £ 4 (t) with 
t = 2.5n. 

We focus first on the range 3.8 < r < 3.94 with increments of 0.00025. The "exact" 
maximal Lyapunov exponent calculated by solving the analytic linearised flow for a 
noise free trajectory is shown in Figure Efa) and serves again as a benchmark to 
compare maximal Lyapunov exponents with our test. 

In Figure Efb,c) we show the results for noise-free data for the direct method and 
for our method. The methods appear to work equally well here. It is interesting 
to note that at r = 3.8175 both methods indicate regular behaviour whereas the 
"exact" maximal Lyapunov exponent is positive at that value of r indicating chaotic 
behaviour. As a matter of fact there is a periodic window close by at r = 3.817525 
which we checked by using the "exact" maximal Lyapunov exponent. 

In FigureEfb), we have taken the same vertical range (from to 0.3) as was natural 
in Figure Efa). Note the poor convergence of the maximal Lyapunov exponent for 
the direct method. 

The results for i.i.d. measurement noise using a noise-level of 10% are shown in 
Figure |H1 The methods seem to perform roughly on an equal level. Whereas our 
method misses the periodic window at r = 3.83025, the direct method misses the 
periodic window at r = 3.88025. 
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Figure 5: Plots of the maximal Lyapunov exponent and K versus r for the 8- 
dimensional Lorenz 96 system (|5.1|) in the regime 3.8 < r < 3.94 using N = 10, 000 
noise-free data points, (a): "Exact" maximal Lyapunov exponent, (b): Maximal 
Lyapunov exponent, direct method [T2]. (c): Our test. 



14 



0.3 



(a) 



(b) 



3.8 3.82 3.84 3.86 3.88 3.9 3.92 3.9 



3.8 3.82 3.84 3.86 3.88 3.9 3.92 3.94 



Figure 6: Plots of the maximal Lyapunov exponent and K versus r for the 8- 
dimensional Lorenz 96 system (|5.1jl in the regime 3.8 < r < 3.94 using N = 
10,000 data points with 10% measurement noise, (a): Maximal Lyapunov, direct 
method 02]. (b): Our test. 
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Next, we consider the quasiperiodic regime 5.25 < r < 5.5 increasing in increments 
of 0.0005. The "exact" maximal Lyapunov exponent is shown in Figure [7(a). In 
Figure |7Jb,c), we show the results for noise-free data for the direct method and for 
our method. Note the difference in the vertical ranges in Figure |7Ja,b). Whereas the 
K test shows a distinction of regular quasiperiodic motion and chaotic motion near 
r = 5.25 the direct method cannot properly resolve this. This becomes more drastic 
in the presence of noise. 

In Figure |H1 we show the results for the quasiperiodic regime when 10% measure- 
ment noise is added. Here the distinction near r = 5.25 is not possible at all for 
the direct method. The K test shows a faint signature close to r = 5.25 but the 
chaotic interval 5.269 < r < 5.27 is clearly visible. Moreover, the direct method has 
a spurious chaotic interval 5.311 < r < 5.312. In addition, the direct method fails 
as an absolute test. For example the periodic region neighbouring the chaotic peak 
at r = 5.4535 is almost indistinguishable from an absolute point of view. The Lya- 
punov exponent of the periodic regime is about 80% the value of the chaotic peak. 
In contrast our method does not have any spurious regular or chaotic intervals. Also 
it works well as an absolute test and the difference between values of K indicating 
regular dynamics and values of K indicating chaotic dynamics are well separated. 
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Figure 7: Plots of the maximal Lyapunov exponent and K versus r for the 8- 
dimensional Lorenz 96 system (|5.1|) in the quasiperiodic regime 5.25 < r < 5.5 using 
N = 10,000 noise-free data points, (a): "Exact" maximal Lyapunov exponent, (b): 
Maximal Lyapunov exponent, direct method (c): Our test. 
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Figure 8: Plots of the maximal Lyapunov exponent and K versus r for the 8- 
dimensional Lorenz 96 system (|5.1jl added measurement noise in the quasiperiodic 
regime 5.25 < r < 5.5 using N = 10, 000 data points with 10% measurement noise, 
(a): Maximal Lyapunov exponent, direct method |12j . (b): Our test. 
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6 Summary 



We have presented a modification of our test for chaos [3] which better handles mod- 
erate amount of contaminated data. In particular we have simplified the diagnostic 
equation, and we now compute the asymptotic growth rate K as the median value 
over realizations over many different values of c. 

We have shown that our test is much better at coping with measurement noise than 
tangent space methods. Moreover, in the case of high dimensional ODE's we found 
that our method compares well also with the direct method proposed by Rosenstein et 
al. ^2] in the noise-free case, and shows much better results in the case of underlying 
quasiperiodic dynamics in the presence of measurement noise. Our test works well as a 
relative test and also as an absolute test for moderately short time series contaminated 
by measurement noise. 

Since our test applies directly to the time-series data, bypassing the need for 
phase space reconstruction, we anticipate that the advantages of our method will 
be magnified for data from partial differential equations and from real experiments. 
(However, the comparison is complicated by the increased difficulty in generating the 
data and the fact that often in these situations the "correct" answer is not known 
beforehand, which is why a reliable test is required in the first place.) This will be 
the subject of future work. 
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